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Abstract. We report on two quantitative, morphological estimators of the filamentary structure of the Cosmic Web, the so- 
called global and local skeletons. The first, based on a global study of the matter density gradient flow, allows us to study 
the connectivity between a density peak and its surroundings, with direct relevance to the anisotropic accretion via cold flows 
on galactic halos. From the second, based on a local constraint equation involving the derivatives of the field, we can derive 
predictions for powerful statistics, such as the differential length and the relative saddle to extrema counts of the Cosmic web 
as a function of density threshold (with application to percolation of structures and connectivity), as well as a theoretical 
framework to study their cosmic evolution through the onset of gravity-induced non-linearities. 

Keywords: cosmology, larges scales structures, topology 
PACS: 95.35.+d, 95.36.+X 98.80.-k, 98.80.-k 98.80.Jk 

Over the course of the last decades, our understanding of the extragalactic universe has undergone a paradigm 
shift: the description of its structures has evolved from being (mostly) isolated to being multiply connected both on 
large scales, cluster scales and galactic scales. This interplay between large and small scales is driven in part by the 
scale invariance of gravity which tends to couple dynamically different processes, but also by a the strong theoretical 
prejudice associated with the so-called concordant cosmological model |T]. This model predicts a certain shape for the 
initial conditions, leading to a hierarchical formation scenario, which predicts the formation of the so-called Cosmic 
Web, the most striking feature of matter distribution on megaparsecs scales in the Universe. This distribution was 
confirmed, observationally, more than twenty years ago by the first CfA catalog and later by the SDSS |3] and 
2dFGRS |4] catalogs. On these scales, the "Cosmic Web" picture (sl] relates the observed clusters of galaxies, and 
filaments that link them, to the geometrical properties of the initial density field that are enhanced but not yet destroyed 
by the still mildly non-linear evolution |@] . The analysis of the connectivity of this filamentary structure is critical to 
map the very large scale distribution of our universe to establish, in particular, the percolation properties of the Web 

a 

On intermediate scales, the paradigm shift is sustained by panchromatic observations of the environment of galaxies 
which illustrate sometimes spectacular merging processes, following the pioneer work of e.g. |0] (motivated by 
theoretical investigations such as [IQD . The importance of anisotropic accretion on cluster and dark matter halo scales 
flTl [T2[ [T3I1 is now believed to play a crucial role in regulating the shape and spectroscopic properties of galaxies. 
Indeed it has been claimed (see e.g. (3, [HI]) that the geometry of the cosmic inflow on a galaxy (its mass, temperature 
and entropy distribution, the connectivity of the local filaments network, etc.) is strongly correlated to its history 
and nature. Specifically, simulations suggest that cold streams violently feed high redshift young galaxies with a 
vast amount of fresh gas, resulting in very efficient star formation. One of the puzzles of galaxy formation involves 
understanding how galactic disks reform after minor and intermediate mergers, a process which is undoubtedly 
controlled by anisotropic gas inflow. 

Recently, |[T6ll presented a method to compute the full hierarchy of the critical subsets of a given density field. This 
approach connects the study of the filamentary structure to the geometrical and topological aspects of the theory of 
gradient flows f[% . In this paper, we focus on the connectivity of the corresponding network. Specifically, since galaxy 
formation seems to be geometrically regulated by the accretion of cold gas from the LSS, how can we use the skeleton 
to characterize this anisotropic accretion, and predict its evolution through perturbation theory ? 

Let us first qualitatively introduce the two operating skeleton extraction algorithms, and present our findings 




FIGURE 1. Left: an illustration of the local condition of the skeleton on a given ridge: on such a critical line, the gradient is 
parallel to the largest eigenvalue of the Hessian (corresponding to the least curvature). This defines a degenerate point process which 
provides means of extending the critical point theory to critical lines. Right: the fully connected skeleton algorithm illustrated in 
2D. Following the panels anti-clockwise, the gradient of a given field defines neighboring regions, which correspond to the set of 
pixels which lead to the same minimum (a so-called void-patch). The edge of these void-patches defines a set of lines which connect 
critical points together: the global skeleton. 




FIGURE 2. Left: the 3D void-patch of a simulation of the cosmological density density field in a 50h~ Mpc box with the 
gadget -2 N-body code. This void-patch segmentation was computed from a 128^ pixels sampling grid smoothed over 5 pixels 
2h~^ Mpc). Right: the corresponding skeleton. The skeleton's different colours represent the index of the peak-patches (i.e., 
the void-patches of minus the density field) which provide by construction the natural segmentation of filaments attached to the 
different clusters. 



regarding the connectivity of random fields. We will then summarize the underlying statistical theory, first for Gaussian 
random fields, then for non-Gaussian fields, which allows us to explain qualitatively the cosmic evolution of the 
connectivity of dark matter halos. In doing so we will also demonstrate how the skeleton applied to the large scale 
structure of the universe could be used to track ^^dm and D{z). 



The skeleton: algorithms 



We may qualitatively define the skeleton as the 3D analog of ridges in a mountainous landscape. It consists of the 
special lines which connect saddles and peaks together. Mathematically these are critical solutions of the gradient flow 
(r = Vp) between critical points. Here we want to construct such 3D ridges to trace the filaments of the large scale 
structures. Two venues have been explored recently: the so called local skeleton Hslll^], which defines a (degenerate) 
point-like process corresponding to the zeros of a set of functions; the idea is that on a ridge the gradient should be 
along the direction of least curvature (see figure [T] left panel), which translates into the set of 3 equations: 

S = (VVp-Vp)x Vp = 0, (1) 

where Vp and VVp are respectively the gradient and the Hessian of the field, p and Ai > A2 > A3 its eigenvalues, with 
the condition Ai + A2 < 0, which ensures that the critical line is a ridge. Alternative conditions on the eigenvalues of the 
Hessian may be applied to pick up the whole set of critical lines, see figure[3] Following |[l9ll . this local skeleton and its 
properties provide an alternative description to classical approaches of galaxy clustering (powerspectrum, bispectrum 
etc.) and attempt to achieve data compression via a mathematical description of the morphology of the cosmic web. 

An alternative, global definition ll^ [l^ is to define it as the border between N void-patches, where N is the 
dimensionality and void-patches are the attraction basins of the (anti-)gradient of the field (see figure [TJ right panel in 
2D, and figure [2l in 3D). The actual implemented algorithm is based on a watershed technique and uses a probability 
propagation scheme to improve the quality of the segmentation (see also ifTsIl ). It can be applied within spaces of 
arbitrary dimensions and geometry. The corresponding recursive segmentation yields the network of the primary 
critical lines of the field: the fully connected skeleton that continuously link maxima and saddle-points of a scalar field 
together. Both constructions are of interest: the local formulation provides means of conducting detailed statistical 
investigation of the corresponding degenerate point-like process, while the global skeleton yields a totally connected 
network of lines. This critical set of lines is a compact description of the geometry of the field, richer than the 
knowledge of the critical points alone. For the purpose of this review, it also allows us to explore the connectivity 
of peaks, both from the point of view of the number of connections to other peaks, and of the number of skeleton 
branches leaving a given maximum. Conversely, the local skeleton formulation allows us to extend the BBKS theory 
1I2TI] of peaks in the context of critical lines to investigate their cosmic evolution, see below. Figure [3] compares the two 
sets of lines in the neighbourhood of a given peak. 



The skeleton: connectivity 

Let us now make use of this algorithm to explore the connectivity of the corresponding network. A set of two 
dimensionnal Gaussian random fields (GRF) is produced, and for each of them, the peak patch and the skeleton was 
generated and its connectivity computed, following the second prescription described in 1 16]: we chose here to smooth 
and label the skeleton having fixed the field extrema, (i.e. not the bifurcation points of the skeleton, see Figured (top 
left panel)). By fixing the extrema of the field, one ensures that the skeleton subsets that link these extrema are treated 
independently. 

Let us first focus on the degree of the peakpatch, defined as the number of saddle points within a given patch which 
connects the skeleton of one patch to its neighbours. Hence the connectivity count reflects the number of maxima 
a given peak is connected to (or in the language of graph theory, the most likely degree of the vertices: fi = (n)). 
Figure m (top right panel) displays this PDF as a function of n; in particular it allow us to compute {n), the mean 
n which is found to be independent of 7, the shape parameter of the underlying GRF. The bottom left panel shows 
the corresponding normalized PDF(^, T]) versus the contrast T]. As expected the number of connections increases with 
contrast. Indeed, near the maximum at high contrast all eigen values tend to become equal ll22ll . Therefore all incoming 
directions become possible. 

Let us now take the other (astrophysically relevant) perspective and count the actual number of filaments incident 
onto a given maximum. This local (intra patch) degree is in fact equal to the number of saddle points minus the 
number of bifurcation within the peakpatch: ^max = ^saddle — ^bifurcation • The bottom right panel of figure IH shows 
the corresponding PDF(/2niax)- Note that this distribution is almost symmetrical, centered at (^max) = 3. Figure [5] (ieft 
panel) presents the same distribution in 3D, which is very skewed and presents a sharp mode at 3. Its cosmic evolution, 
measured in CDM dark matter simulations with and without dark energy is qualitatively shown on the right panel. Our 
purpose in the rest of this paper is to explain qualitatively this cosmic trend by 1) deriving the statistics of bifurcation 
points within each patch and 2) predicting the cosmic evolution of saddle points and peaks. 




FIGURE 3. An example of a generic patch of a 2D field. The underlying isocontours correspond to the density field. The thin 
gold lines show the gradient lines of the field. The blue lines represent the local set of critical lines, given by the solution of equation 
IT]). Primary lines are shown in solid and secondary lines are dashed. The green lines correspond to global critical lines: the skeleton 
and the anti- skeleton, which delineate a special bundle of gradient lines at resp. the intersection of a peak-patch and a void-patch 
CI . The primary local lines follow fairly well the gradient lines, noticeably near the saddle points, where the "stiff" approximation 
holds best. In contrast, the approximation worsens for the secondary critical lines. The main distinction between the global and 
local skeletons is that the global one follows everywhere the smooth gradient line that uniquely connects a maximum to a saddle, at 
the cost of deviating from being exactly on the ridge (see how in the vicinity of the minimum at the bottom, the right green line does 
not follow the valley). The local skeleton tries to delineate the ridges as far from extrema as possible, but then the lines that follow 
this local procedure from different extrema do not meet and have to rather suddenly reconnect. A particularly striking example is 
the loop on the right hand side. Bottom left medallion: the neighbourhood of a local critical line (thick blue line). Thin lines are 
isocontours of the field. Three sample points are investigated in detail. The signature, orientation and the magnitude of the local 
Hessian are represented by the golden shapes. Near the maximum on the right edge, the signature of the eigenvalues of the Hessian 
is (-1,-1), which is shown by ellipses oriented according to eigen-directions with longer semi-axis along the direction of the least 
curvature. At the leftmost point the eigenvalue signature is "saddle-like", (1,-1), which is represented by a pair of hyperbolae, also 
oriented with respect to eigen-directions. By definition, on the critical line the gradient of the field Vp, shown by red arrows, is 
aligned with one of the eigen-directions (i.e the axis of the ellipse or hyperbola in the graph). The light cyan arrows are the tangent 
vectors to the critical line oc £ . V5, while stiff approximation to them would be parallel to the gradient. The direction of the critical 
line is close to the gradient when it follows the ridge near the maximum, but slides at an angle in the "saddle-like" region, before 
joining the corresponding saddle extremal point. Note that the gradient line that takes us to the same saddle as a segment of the 
global skeleton (dashed green) does not follow the ridge too closely in this instance. 



The skeleton: statistics 



As illustrated first in 3D on Figure[7l the critical condition given by equation ([T]) defines 3 isosurfaces corresponding 
to its components, and the local skeleton direction is simply given by the cross products of two of its normals. Hence 
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FIGURE 4. Top left panel: A typical (gray coded) 2D peak-patch of a Gaussian random field together with the degree of each 
vertex of the skeleton graph; for instance, within the red star the mean degree is 4, while the corrected degree is 3 within the blue 
star, as one bifurcation point is present; top right panel: the corresponding PDF of this degree, for different values of the shape 
parameter 7 as labeled; the mean degree is 4 and corresponds to the theoretical expectation of 2;isaddie/'^max; bottom left panel: 
the (log) differential PDF of the degree as a function of contrast: as expected, the denser peaks are more connected; interestingly, 
the distribution seems Gaussian at fixed n even though its marginal isn't; bottom right panel: the corrected PDF of the number of 
branches connecting onto a given peak, which involves substracting the number of bifurcation points within the patch to the number 
of saddle points on the edge of the patch. 



the differential length (per unit volume) is simply given by the statistical expectation 

(f'Xkd^Xkld^^Xklm X Vsj\P{r{,Xk,Xkl,Xklm)8D {s\Xk,Xkl,Xklm)) & {s^ {Xk.Xkl.Xklm)) , (2) 

where ergodicity allowed us to replace volume average by ensemble average over the statistical distribution, 
P{T\^Xk-,Xki-,Xkim)^ of the successive (a-dimensional) derivatives, G\Xk = Vy^p, 02x^1 = ^k^ip^ <^?>^klm = ^mSi^kP of 
the field, p. The scale factor 1 /R^^ = of /erf in equation Q follows dimensionally (given oj ^jV^ = VS\ GqT] = p). 
Here the variances, O"^, obey Gp = In^^'^ /T[N /2] Jq k^PP{k)k^~^dk in N dimensions, (see ||2Q|i for a formal derivation 
of equation ©). 
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FIGURE 5. Left panel: the 3D corrected PDF of the number of branches connecting onto a given peak for a GRF: the mode is 
at 3; insert: the typical neighbourhood entropy map of a massive galaxy in the Marenostrum simulation from fisll : right panel: the 
qualitative cosmic evolution of the mean degree in dark matter simulation with and without dark energy; as expected this degree 
decreases with the expansion factor as connectors get washed out by gravitationnal clustering; the rate of decrease seems sensitive 
to the moment when dark energy kicks in for ACDM models (around a ~ 0.5). Note that the absolute amplitude of this mean degree 
is not shown here as we expect the global skeleton to over estimate the mean connectivity in 3D (see main text). 



More generally, for ND critical lines, the N-l independent functions that define the critical condition ([U acquire 
the following reduced form in the eigenframe of the Hessian of the field: = xiXi (Ai — A/) = 0. The measurements 
in 1I23I] found that, over the range of spectral indexes relevant to cosmology, the third derivatives of the field x^im 
play a subdominant nature. In the so called "stiff" approximation we therefore omit the third derivative, effectively 
assuming that the Hessian can be treated as constant during the evaluation of V^. This picture corresponds to 
a skeleton connecting extrema with relatively straight segments. Let us first assume that the underlying field is 
Gaussian. In the stiff approximation, the gradients s\ = Vs^ have just two non-zero components, ^-'i = — A/) 

(which vanishes on the critical line) and s^i = xi — A/), which shows that in this approximation we equate the 
direction of the line with the gradient of the field. Substituting this expression into equation Q and integrating over 
5d(^0 = ^{^i)/ (-^1 (^1 ~ ^i)) we obtain a simple expression for the differential length of the ND-critical lines: 
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where the shape parameter 7 = 0^/(0002) describes the correlation between the field and its second derivatives, Qy is 
a quadratic form in A, and rj which functional form is 



Qr{ri,{^i}) = ri' 
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Equation ^ simply states that the stiff differential length is the expectation of the product of the "smaller" 
eigenvalues, which is quite reminiscent of the classical extremal count result (the expectation of the product of all 
eigenvalues). It also qualitatively makes sense, as the larger the curvarture orthogonal to the skeleton, the more skeleton 
segment one may pack per unit volume. Since the argument of Qy is extremal as a function of 77 when yrj ~ J^, A,, 
the largest contribution at larg e yr] in the integral should arise when A; 0= yr] since near the maximum at high contrast 
all eigen values are equal 1I22I1 . Hence given that Yli<j{^ — Ay) is the measure, the only remaining contribution in the 
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FIGURE 6. The relative differential length, d^/dr] /PDF, measured in simulation of 2D (left) and 3D (right) Gaussian random 
fields with scale invariant power-law spectra versus predictions of the local theory in stiff approximation (solid curves). The spectral 
parameters are 7= 0.71,0.59,0.39 for the 2D and 7= 0.77,0.70,0.60 for the 3D simulations. 



integrand comes from \ Y[i>i ^i\ {^V) ' the dominant term at large T] is given by 
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Note that in 2D the differential length, equation ([3]) can be reduced to a particularly simple closed form 
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and for the integrated skeleton length 
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In other words, one expects to find one segment of skeleton per linear section of ~ (4.27?*). Similarly, in 3D 
one segment of skeleton is found per surface section of ^ (4.657?^). The match of the detailed PDF with the 
corresponding measurements are shown in Figure [6] in 2 and 3D. From the point of view of cosmology, the skeleton 
is invariant w.r.t. any monotonic bias, and traces the denser regions of the field. As shown in equations ^ and ^ its 
statistical description yields a measure of both and Rq, hence on the shape of the underlying powerspectrum on 
the corresponding smoothing scale over which the skeleton was computed. An implementation of this estimate on the 
SDSS catalog was carried by ll24l] and provided constraints on ^^m- 



Bifurcation counts 

In 1I25I] we conjectured that critical lines experience a qualitative change in behaviour in the vicinity of the points 
where either a Hessian eigenvalue orthogonal to the gradient direction vanish, or becomes equal to the one along 
the gradient. The first type corresponds to points where the curvature transverse to the direction of the critical line 
vanishes along at least one axis: typically, in 2D, they mark regions where a crest becomes a trough, or vanish into 
a plateau. The second type correspond to points where the critical lines would split, even though the field does not 
go through an extremum: a bifurcation of the lines occurs along the slope; the occasional skier or mountaineer will 
be familiar with a crest line splitting in two, even though the gradient of the field has not vanished. Namely, if, for 
definiteness, Vp is taken to be along the first eigen-direction, A2 = 0, or A2 = Ai. We called the first case the "sloping 
plateau" as it designates the entering of a flat region, and the second, the "bifurcation" as it designates the places 
of possible reconnection of critical lines. Remarkably, these special points on the critical lines are recovered by the 



FIGURE 7. left panel: the local direction of the skeleton is along the cross product, V^Si x VS2 of the normals, V^Si and VS2, 
of the iso-surfaces, and ^2, defined by the critical condition ([T]). Hence the differential length per unit volume of the skeleton 
is given by expectation of the modulus of this cross product, which can be expressed locally with the JPDF of the field and its 
successive derivatives, right panel: the three types of singular point (V^^^S = 0) on the critical lines (in solid blue: primary; dashed: 
secondary; green: gradient lines of the global skeleton): an extremum (open circle), a "bifurcation" (white square), and "slopping 
plateaux" (black squares). The red star corresponds to what would naively correspond to the bifurcation point. 



formal singular condition zero-tangent vector if VkS^ is evaluated in the stiff approximation. Along the ND critical line 
defined by X2 = . . . = x^v = 0, the nul tangent vector condition, x^~^ Y[i>i h {^\ — h) = gives rise to three classes 
of situations: (i) xi =0 corresponding to extremal points; (ii) one of A/ = corresponding to slopping flattened tubes; 
and (iii) one of A; = Ai, corresponding to an isotropic bifurcation. 

To be more specific, in 2D, the skeleton's singular points correspond to points where Sk = = 0. The number 
density, ^b(^) of singular points below the threshold T] is equal to 
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The gradient of S, evaluated in the stiff approximation, in the Hessian eigenframe has the components — 
X2A1 (Ai — A2) , and S2^^ = xi A2 (Ai — A2) , and involves only second derivatives of the field. Let us consider the critical 
line that corresponds to the X2 = condition in the Hessian eigenframe. Then sf^^ vanishes everywhere along this line. 
The requirement ^2^^ = has a solution at the extremal points, xi = 0, but also in two other cases, namely A2 = or 
A2 = Ai, that we conjectured to be of interest. The second situation (isotropic Hessian) has a number density given by 
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We note that the number density of "bifurcation" points is proportional just to the PDF of the field and, consequently 
the bifurcation points are as frequent in the regions of high field values as in the low ones. Finally, in 3D, we expect 
the finite resolution bifurcation branch (connecting maxima to bifurcation points) to become bifurcation plates; in 
turn the boundary of these plates will be counted as critical lines by our voidpatch algorithm; we therefore expect 
that the algorithm numerically over estimate the degree of peakpatches. Note nonetheless that the 3D counterpart of 
equation © should allow us to correct for the number of bifurcations within each patch and compute the statistics of 
the number of branches connected to a given maximum. 



Departure from Gaussianity 



While the Gaussian limit provides the fundamental starting point in the study of random fields ll26[[27[ [2lll. non- 
Gaussian features of the CMB and the large scale structure (LSS) fields are also of great interest. CMB inherits high 




FIGURE 8. The number of extrema (left) and the Euler characteristic (right) for 3D gravitational collapse. In dimensional 
units both quantities are given per rI volume. The shaded bands correspond to 2(7 variations over the mean measurements of 
25 realizations, while the curves give the 3^^ and 4* order predictions. The Gaussian prediction is shown as a thin line. 



level of Gaussianity from the initial fluctuations, and small non-Gaussian deviations may provide a unique window into 
the details of processes in the early Universe. The gravitational instability that nonlinearly maps the initial Gaussian 
inhomogeneities in matter density into the LSS, on the other hand, induces strong non-Gaussian features culminating 
in the formation of collapsed, self-gravitating objects such as galaxies and clusters of galaxies. At supercluster scales 
where non-linearity is mild, the non-Gaussianity of the matter density is also mild, but still essential for quantitative 
understanding of the filamentary Cosmic Web fy\ in-between the galaxy clusters. 

In order to extend the result of the previous section, following [28], let us develop the equivalent of the Edgeworth 
expansion for the JPDE of the field variables that are invariant under coordinate rotation. Such distribution can be 
obtained directly from general principles: the moment expansion of the non-Gaussian JPDE corresponds to the 
expansion in the set of polynomials which are orthogonal with respect to the weight provided by the JPDE in the 
Gaussian limit. Thus, the problem is reduced to finding such polynomials for a suitable set of invariant variables. The 
rotational invariants that are present in the problem are: the field value T] itself, the modulus of its gradient, = Y^i^i^ 
and the invariants of the matrix of the second derivatives Xij. A rank N symmetric matrix has N invariants with respect 
to rotations. The eigenvalues provide one such representation of invariants, however they are complex algebraic 
functions of the matrix components. An alternative, more useful representation is given by the linear combination of 
the polynomial invariant. Is (where the linear invariant is the trace, Ii = the quadratic one is h = lli<j K^j and the 

N-th order invariant is the determinant of the matrix. 1^ = Ui ^i) J\=h, Js>2 = If - Ep=2 ( -N^Cf /i^-^) /^P\~^h ' 
where Js>2 are (renormalized) coefficients of the characteristic equation of the traceless part of the Hessian and are 
independent in the Gaussian limit on the trace J\ . Let us consider again here the 2D case explicitly. Introducing 
= (t] + 7/1)/ — in place of the field value T] we find that the 2D Gaussian JPDE G2v>(J^^q^ ^h^^i)^ normalized 
over d^dc^dJ\dJ2, has a fully factorized form in these variables Used as a kernel for the polynomial expansion, G2D 
leads to a non-Gaussian rotation invariant JPDE in the form of the direct series in the products of Hermite (Hi), for ^ 
and /i, and Laguerre (LO, for (f^ and J 2, polynomials: 
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where Y!i^j^u=o^^~^ stands for summation over all combinations of non-negative /, j, I such that i-\-2j -\-k-\-2l adds 
to the order of the expansion term n. A similar expression holds in 3D. The coefficients of the expansion are "centered" 
moments, given by the differences of the actual moments and their Gaussian limit (the latter vanishing when either / 
or k is odd): {^q^^ Ji^J2^)c = {Cq^^Ji^J2^) - {i - l)\\{k- l)\\jU\ . Once the JPDE is known it is straightforward to 
compute any expectation of the field which involve algebraic combinations of the invariants such as the differential 



length, d^^^ /drf (equation ©), extrema counts, or the Euler characteristic. For instance, the Euler characteristic 
can be computed completely (since there are no sign constraints on the eigenvalues of the Hessians) by noting that 
In depends only linearly on Js>2 (e.g., I3 = {Ji^ — 3/1/2 + 2/3) /27 in 3D), hence all terms in JPDF of higher order 
in Js>2 do not contribute. The 2D and 3D results can be combined in a very compact form if one re-expresses the 
"centered" moments back in terms of the field T] itself and the invariants Is 
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where i = 0,s = N terms have been combined into the boundary term oc ^(—00) fixed by the topology of the manifold, 
and should be omitted from the further sum. Now if the departure from Gaussianity is induced by gravitationnal 
clustering, the cumulants occuring in equations ([TOl) and (fTTI) can be computed in the context of perturbation theory 
[29, 30], will scale like the growth factor D{z), and can be used to constrain the dark energy equation of state via 3D 
galactic surveys, or shed light on the physics of the early Universe through 2D CMB maps. Regarding the connectivity, 
equation (fTTI) is illustrated on figure [8] in this context. In particular, the non-linear evolution of the number of saddle 
points and the number of peaks is accurately predicted, which in turn should allow us to predict the mean degree, (n) 
of the peaks within each peak-patch given that each saddle-point connects to two peaks: {n{z)) = 2^saddie(^)/^max(z). 
When a non-Gaussian extension of equation ^ is derived we should also be in a position to predict the non-linear 
evolution of the number of connecting streams on dark matter halos. 
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